{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quadratic Programming\n",
    "cvxopt 中的二次规划标准式：\n",
    "$$\n",
    "\\begin{align}\n",
    "\\min_{x}\\ &\\frac{1}{2}x^TPx+q^Tx \\\\\n",
    "subject\\ to\\ &Gx\\preceq h \\\\\n",
    "         &Ax=b\n",
    "\\end{align} \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cvxopt as co\n",
    "co.solvers.options['show_progress'] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_hard_margin_SVM(X, y):\n",
    "    \"\"\"\n",
    "    Linear Hard Margin SVM by QP\n",
    "    Args:\n",
    "        X: 数据\n",
    "        y: 标签\n",
    "    \n",
    "    Returns:\n",
    "        b: intercept\n",
    "        w: 特征权重\n",
    "    \"\"\"\n",
    "    P = np.block([\n",
    "        [0.0,                                    np.zeros_like(X[0])],\n",
    "        [np.zeros_like(X[0].reshape((-1,1))),    np.eye(len(X[0]))  ]\n",
    "    ])\n",
    "    q = np.zeros((len(X[0])+1,1), dtype=np.float64)\n",
    "    G = -y.reshape((-1,1)) * np.hstack((np.ones((len(y), 1), dtype=np.float64), X))\n",
    "    h = -1*np.ones_like(y, dtype=np.float64).reshape((-1,1))\n",
    "    \n",
    "    result =  np.array(co.solvers.qp(co.matrix(P), co.matrix(q), co.matrix(G), co.matrix(h))['x'])\n",
    "    return result[0], result[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0]\n",
      " [2 2]\n",
      " [2 0]\n",
      " [3 0]] [-1 -1  1  1]\n",
      "[[ 0  0]\n",
      " [ 2  2]\n",
      " [-2  0]\n",
      " [-3  0]]\n"
     ]
    }
   ],
   "source": [
    "X = np.array([[0,2,2,3], [0,2,0,0]]).T\n",
    "y = np.array([-1,-1,1,1])\n",
    "print(X, y)\n",
    "print(-y.reshape((-1,1))*X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b = [-1.00000001], \n",
      "w = [[ 1.00000001]\n",
      " [-1.00000001]]\n"
     ]
    }
   ],
   "source": [
    "b, w = linear_hard_margin_SVM(X, y)\n",
    "print('b = {}, \\nw = {}'.format(b, w))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_hard_margin_SVM_dual(X, y):\n",
    "    \"\"\"\n",
    "    Linear Hard Margin SVM Dual by QP\n",
    "    Args:\n",
    "        X: 数据\n",
    "        y: 标签\n",
    "    \n",
    "    Returns:\n",
    "        alpha: 拉格朗日乘子系数\n",
    "        sv_index: support vector index\n",
    "        b: intercept\n",
    "        w: 特征权重\n",
    "    \"\"\"\n",
    "    P = y.reshape((-1,1))*(X.dot(X.T))*y.reshape((1,-1)).astype(np.float64)\n",
    "    q = np.ones_like(y) * -1.0\n",
    "    G = -1*np.eye(len(y),dtype=np.float64)\n",
    "    h = np.zeros_like(y, dtype=np.float64).reshape((-1,1))\n",
    "    A = y.reshape((1,-1)).astype(np.float64)\n",
    "    b = 0.0\n",
    "    \n",
    "    alpha =  np.array(co.solvers.qp(co.matrix(P), co.matrix(q), co.matrix(G), co.matrix(h), co.matrix(A),co.matrix(b))['x']).reshape((-1,))\n",
    "    sv_index = np.where(alpha>1e-6)[0]\n",
    "    \n",
    "    w = (alpha[sv_index]*y[sv_index]).dot(X[sv_index])\n",
    "    n = np.random.choice(sv_index)\n",
    "    b = y[n] - w.dot(X[n])\n",
    "    return alpha, sv_index, w, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha:  [5.00000000e-01 5.00000006e-01 9.99999998e-01 7.80987988e-09]\n",
      "w:  [ 0.99999999 -1.00000001]\n",
      "b:  -0.9999999478691657\n"
     ]
    }
   ],
   "source": [
    "alpha, _, w, b = linear_hard_margin_SVM_dual(X, y)\n",
    "print('alpha: ', alpha)\n",
    "print('w: ', w)\n",
    "print('b: ', b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha:  [2.49999916e-01 8.17639192e-08 2.49999998e-01]\n",
      "w:  [0.49999975 0.49999975]\n",
      "b:  -1.9999985061835712\n"
     ]
    }
   ],
   "source": [
    "X, y = np.array([[3,3],[4,3],[1,1]]), np.array([1,1,-1])   # 李航老师书的例子\n",
    "alpha, _, w, b = linear_hard_margin_SVM_dual(X, y)\n",
    "print('alpha: ', alpha)\n",
    "print('w: ', w)\n",
    "print('b: ', b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
